
close all;
clear all;
fclose all;

%create data structure
HCR=struct('CC',[],'MLI2',[],'Globular',[],'Lugaro',[],'Golgi',[]);

% path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\4. position measurements\';
path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\OTRcrex Ai14 Slc6a5 Aldh1a3 12.7.20\Processing\4. position measurements\';

dir_to_search=[path1];


% get the folder contents
d = dir(path1);
% remove all files (isdir property is 0)
dfolders = d([d(:).isdir]==1);
% remove '.' and '..'
dfolders = dfolders(~ismember({dfolders(:).name},{'.','..'}));

for SliceFolder=1:numel(dfolders) %slices
    temp_dir=[dfolders(SliceFolder).folder '\' dfolders(SliceFolder).name '\'];
    formatpattern = fullfile(temp_dir, '/*.mat');
    
    dInfo=dir(formatpattern);
    
    
    GolgiloadFileIndex=nan;
    CCloadFileIndex=nan;
    GlobularloadFileIndex=nan;
    MLI2loadFileIndex=nan;
    LugaroloadFileIndex=nan;
    
for j=1:numel(dInfo) %Cell types
    %make logical for CC,globular etc in dInfo to load the correct file
    if ~isempty(strfind(dInfo(j).name,'CC'))
         CCloadFileIndex=j;
    else
         if   ~isempty(strfind(dInfo(j).name,'Globular'))
         GlobularloadFileIndex=j;
         else
             if ~isempty(strfind(dInfo(j).name,'MLI2'))
                  MLI2loadFileIndex=j;
             else
                 if ~isempty(strfind(dInfo(j).name,'Golgi'))
                    GolgiloadFileIndex=j;
                 else
                     if ~isempty(strfind(dInfo(j).name,'Lugaro'))
                        LugaroloadFileIndex=j;
                     end
                 end
             end
         end
    end
end
      
%% Load CCs
if ~isnan(CCloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
    
    load([dInfo(CCloadFileIndex).folder '\' dInfo(CCloadFileIndex).name]);
    HCR(SliceFolder).slice=CellsName(1:4);
    
    if exist('PC')
        HCR(SliceFolder).PC=PC;
        HCR(SliceFolder).ML=ML;
    end
    
    %calculate lobule length using PC layer boundaries and PC this is in
    %pixels of half resolution files of HCR from slidescanner
    if lobulePCboundaries(2)<lobulePCboundaries(1); % outline was traced in wrong direction, lobule X first.
        
        temp=[size(PC,1) lobulePCboundaries 1];
        temp=fliplr(temp); %need to be incremental for indexing purposes, will reverse again later after i'm done calculating length.
    else
        temp=[1 lobulePCboundaries size(PC,1)];
    end
    %
    
    for k=1:8
        lobuleSegments=PC([temp(k):temp(k+1)],:);
        d = diff(lobuleSegments);
        total_length(k) = sum(sqrt(sum(d.*d,2)));
    end
    
    if lobulePCboundaries(2)<lobulePCboundaries(1) % outline was traced in wrong direction, lobule X first reverse it again to get correct order
        total_length=fliplr(total_length);
    end
    
    HCR(SliceFolder).lobuleLength=total_length; % lobule I, II, IV-V,VI,VII,VII,IX,X
    
    % now onto CC stuff
    %asign lobule to specific cells
    lobuleID=[1 2 5 6 7 8 9 10];
    lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)];
    Cellsperlobule=diff(lobuleCellIdx);
    HCR(SliceFolder).CCsperlobule=Cellsperlobule;
    
    clear lobule
    for l=1:8
        
        
      
        lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l);
        end
    end
    
 
    
    
 
    CCdata=[Cells distancePCL distanceML MLlength lobule']
    HCR(SliceFolder).CC=CCdata;
    clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx CCdata 
end
%% Load Globular
if ~isnan(GlobularloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
     load([dInfo(GlobularloadFileIndex).folder '\' dInfo(GlobularloadFileIndex).name]);
    lobuleID=[1 2 5 6 7 8 9 10];
    lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)]; %complete missing values (start and end).
    Cellsperlobule=diff(lobuleCellIdx); % calculate cells per lobule in order of lobuleID.
    HCR(SliceFolder).Globularperlobule=Cellsperlobule;
    clear lobule
    for l=1:8
        lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
    end
    Globulardata=[Cells distancePCL distanceML MLlength lobule'];
    HCR(SliceFolder).Globular=Globulardata;
    clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx Globulardata
end
%% load Golgi
if ~isnan(GolgiloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
    load([dInfo(GolgiloadFileIndex).folder '\' dInfo(GolgiloadFileIndex).name]);
    lobuleID=[1 2 5 6 7 8 9 10];
    lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)]; %complete missing values (start and end).
    Cellsperlobule=diff(lobuleCellIdx); % calculate cells per lobule in order of lobuleID.
    HCR(SliceFolder).Golgiperlobule=Cellsperlobule;
    clear lobule
    for l=1:8
        lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
    end
    Golgidata=[Cells distancePCL distanceML MLlength lobule'];
    HCR(SliceFolder).Golgi=Golgidata;
    clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx Golgidata Cellsperlobule
end
%% load MLI2
if ~isnan(MLI2loadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
    load([dInfo(MLI2loadFileIndex).folder '\' dInfo(MLI2loadFileIndex).name]);

    lobuleID=[1 2 5 6 7 8 9 10];
    lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)]; %complete missing values (start and end).
    Cellsperlobule=diff(lobuleCellIdx); % calculate cells per lobule in order of lobuleID.
    HCR(SliceFolder).MLI2perlobule=Cellsperlobule;
    for l=1:8
        lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
    end
    MLI2data=[Cells distancePCL distanceML MLlength lobule'];
    HCR(SliceFolder).MLI2=MLI2data;
    clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx MLI2data
end
%% load lugaro
if ~isnan(LugaroloadFileIndex) %only run the section if there is a mat file for the corresponding Cell type, otherwise skip to next section
    load([dInfo(LugaroloadFileIndex).folder '\' dInfo(LugaroloadFileIndex).name]);
    lobuleID=[1 2 5 6 7 8 9 10];
    lobuleCellIdx=[1 lobuleCellIdx size(Cells,1)]; %complete missing values (start and end).
    Cellsperlobule=diff(lobuleCellIdx); % calculate cells per lobule in order of lobuleID.
    HCR(SliceFolder).Lugaroperlobule=Cellsperlobule;
    clear lobule
    for l=1:8
        lobule(lobuleCellIdx(l):lobuleCellIdx(l+1))=lobuleID(l)
    end
    Lugarodata=[Cells distancePCL distanceML MLlength lobule'];
    HCR(SliceFolder).Lugaro=Lugarodata;
    clear Cells distancePCL distanceML MLlength lobule lobuleCellIdx Lugarodata
end





%% Plot shit
hold on
CCPCLdistance=[];
CCMLlength=[];
GlobularPCLdistance=[];
GlobularMLlength=[];
MLI2PCLdistance=[];
MLI2MLlength=[];
GolgiPCLdistance=[];
GolgiMLlength=[];
LugaroPCLdistance=[]
LugaroMLlength=[]

 GlobularSliceCount=0;
 GolgiSliceCount=0;
 MLI2SliceCount=0;
 CCSliceCount=0;
 LugaroSliceCount=0;


for SliceFolder=1:numel(HCR)
    if  ~isempty(HCR(SliceFolder).CC)
        CCPCLdistance=[CCPCLdistance; HCR(SliceFolder).CC(:,3)];
        CCMLlength=[CCMLlength ; HCR(SliceFolder).CC(:,5)];
        CCSliceCount=CCSliceCount+1
    end
    
    if  ~isempty(HCR(SliceFolder).Globular)
        GlobularPCLdistance=[GlobularPCLdistance; HCR(SliceFolder).Globular(:,3)];
        GlobularMLlength=[GlobularMLlength ; HCR(SliceFolder).Globular(:,5)];
        GlobularSliceCount=GlobularSliceCount+1
    end
    
    if  ~isempty(HCR(SliceFolder).Golgi)
        GolgiPCLdistance=[GolgiPCLdistance; HCR(SliceFolder).Golgi(:,3)];
        GolgiMLlength=[GolgiMLlength ; HCR(SliceFolder).Golgi(:,5)];
         GolgiSliceCount=GolgiSliceCount+1
    end
    
    if  ~isempty(HCR(SliceFolder).MLI2)
        MLI2PCLdistance=[MLI2PCLdistance; HCR(SliceFolder).MLI2(:,3)];
        MLI2MLlength=[MLI2MLlength ; HCR(SliceFolder).MLI2(:,5)];
         MLI2SliceCount= MLI2SliceCount+1
    end
    
    if  ~isempty(HCR(SliceFolder).Lugaro)
        LugaroPCLdistance=[LugaroPCLdistance; HCR(SliceFolder).Lugaro(:,3)];
        LugaroMLlength=[LugaroMLlength ; HCR(SliceFolder).Lugaro(:,5)];
        LugaroSliceCount=LugaroSliceCount+1;
    end
end

binWidth=0.03
figure;CChist=histogram(CCPCLdistance./CCMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
hold on
Globularhist=histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
Golgihist=histogram(GolgiPCLdistance./GolgiMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
MLI2hist=histogram(MLI2PCLdistance./MLI2MLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
Lugarohist=histogram(LugaroPCLdistance./LugaroMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')


% histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',0.05,'DisplayStyle','stairs', 'EdgeColor', [0 0 0])
% histogram(GolgiPCLdistance./GolgiMLlength,'binwidth',0.05,'FaceAlpha',0.5,'LineStyle', 'none')
% histogram(MLI2PCLdistance./MLI2MLlength,'binwidth',0.05,'FaceAlpha',0.5,'LineStyle', 'none')



%%
color5=[27 192 218]./255;
color1=[0.29, 0.72, 0.35];
color2=[233 83 30]./255; % MLI2 color
 color6=[30, 180, 233]./255; % CCs color from illustrator
color7=[74 184 89]./255; %Globular color
color3=[0 0 0];
color4=[78 233 30]./255;



figure
hold on




GlobularStair=stairs(Globularhist.BinEdges(1:end-1), Globularhist.Values./CCSliceCount)
% Get the (y,x) coordinates from the original inputs or from  
% the xdata and ydata properties of the stair object which will
% always be row vectors.  
bottom = 0; %identify bottom; or use min(b.YData)
x = [GlobularStair.XData(1),repelem(GlobularStair.XData(2:end),2)];
y = [repelem(GlobularStair.YData(1:end-1),2),GlobularStair.YData(end)];
% plot(x,y,'y:') %should match stair plot
% Fill area
area(x,y,'FaceColor',color1,'Facealpha', 0.5,'EdgeColor',color1, 'EdgeAlpha',1,'LineWidth',1)


MLI2Stair=stairs(MLI2hist.BinEdges(1:end-1), MLI2hist.Values./MLI2SliceCount)
% Get the (y,x) coordinates from the original inputs or from  
% the xdata and ydata properties of the stair object which will
% always be row vectors.  
bottom = 0; %identify bottom; or use min(b.YData)
x = [MLI2Stair.XData(1),repelem(MLI2Stair.XData(2:end),2)];
y = [repelem(MLI2Stair.YData(1:end-1),2),MLI2Stair.YData(end)];
% plot(x,y,'y:') %should match stair plot
% Fill area
area(x,y,'FaceColor',color2,'Facealpha', 0.05,'EdgeColor',color2, 'EdgeAlpha',0.7,'LineWidth',1.5)


GolgiStair=stairs(Golgihist.BinEdges(1:end-1), Golgihist.Values./GolgiSliceCount)
% Get the (y,x) coordinates from the original inputs or from  
% the xdata and ydata properties of the stair object which will
% always be row vectors.  
bottom = 0; %identify bottom; or use min(b.YData)
x = [GolgiStair.XData(1),repelem(GolgiStair.XData(2:end),2)];
y = [repelem(GolgiStair.YData(1:end-1),2),GolgiStair.YData(end)];
% plot(x,y,'y:') %should match stair plot
% Fill area
area(x,y,'FaceColor',color3,'Facealpha', 0.05,'EdgeColor',color3, 'EdgeAlpha',0.7,'LineWidth',1.5)

CCStair=stairs(CChist.BinEdges(1:end-1), CChist.Values./CCSliceCount)
% Get the (y,x) coordinates from the original inputs or from  
% the xdata and ydata properties of the stair object which will
% always be row vectors.  
bottom = 0; %identify bottom; or use min(b.YData)
x = [CCStair.XData(1),repelem(CCStair.XData(2:end),2)];
y = [repelem(CCStair.YData(1:end-1),2),CCStair.YData(end)];
% plot(x,y,'y:') %should match stair plot
% Fill area
area(x,y,'FaceColor',color6,'Facealpha', 0.5,'EdgeColor',color6, 'EdgeAlpha',1,'LineWidth',2)


LugaroStair=stairs(Lugarohist.BinEdges(1:end-1), Lugarohist.Values./LugaroSliceCount)
% Get the (y,x) coordinates from the original inputs or from  
% the xdata and ydata properties of the stair object which will
% always be row vectors.  
bottom = 0; %identify bottom; or use min(b.YData)
x = [LugaroStair.XData(1),repelem(LugaroStair.XData(2:end),2)];
y = [repelem(LugaroStair.YData(1:end-1),2),LugaroStair.YData(end)];
% plot(x,y,'y:') %should match stair plot
% Fill area
area(x,y,'FaceColor',[127 63 152]./255,'Facealpha', 0.5,'EdgeColor',[127 63 152]./255, 'EdgeAlpha',1,'LineWidth',2)


vline(0,'--k')
xlim([-1 1])
ylabel('count')

%% Make bar graph of density distribution

%make matrix for the numbers of each cell type (2d Vector: celltype(lobule, slice)
CC_lobules=nan(8,numel(HCR));
Globular_lobules=nan(8,numel(HCR));
Golgi_lobules=nan(8,numel(HCR));
MLI2_lobules=nan(8,numel(HCR));
Lugaro_lobules=nan(8,numel(HCR));

for SliceFolder=1:numel(HCR)
    CC_lobules(:,SliceFolder)=HCR(SliceFolder).CCsperlobule;
    Globular_lobules(:,SliceFolder)=HCR(SliceFolder).Globularperlobule; %Always there
    
    if ~isempty(HCR(SliceFolder).MLI2perlobule)
        MLI2_lobules(:,SliceFolder)=HCR(SliceFolder).MLI2perlobule;
    end
    if ~isempty(HCR(SliceFolder).Golgiperlobule)
        Golgi_lobules(:,SliceFolder)=HCR(SliceFolder).Golgiperlobule;
    end
    if ~isempty(HCR(SliceFolder).Lugaroperlobule)
    Lugaro_lobules(:,SliceFolder)=HCR(SliceFolder).Lugaroperlobule;
    end
    lobule_length(:,SliceFolder)=HCR(SliceFolder).lobuleLength;
end

%% Do bar graphs first without any normalization
figure;
bar(nanmean(CC_lobules,2))
hold on;
for SliceFolder=1:numel(HCR)
scatter([1:size(CC_lobules,1)] ,CC_lobules(:,SliceFolder),'filled')
end
box off

figure;
bar(nanmean(Globular_lobules,2))
hold on;
for SliceFolder=1:numel(HCR)
scatter([1:size(Globular_lobules,1)] ,Globular_lobules(:,SliceFolder),'filled')
end
box off

figure;
bar(nanmean(MLI2_lobules,2))
hold on;
for SliceFolder=1:numel(HCR)
scatter([1:size(MLI2_lobules,1)] ,MLI2_lobules(:,SliceFolder),'filled')
end
box off

figure;
bar(nanmean(Golgi_lobules,2))
hold on;
for SliceFolder=1:numel(HCR)
scatter([1:size(Golgi_lobules,1)] ,Golgi_lobules(:,SliceFolder),'filled')
end
box off

%% make bar plots normalized to PC layer length
%analysis and registration for length of PCL was done in half resolution
%files. Full res scale is 160.67 nm/px, so here it's twice that (321.34)
scale_factor=321.34/1e3; %convert it to microns/px.
lobule_lengthMicron=lobule_length*scale_factor

CC_lobulesNorm=CC_lobules./lobule_lengthMicron*100; %number of cells per 100 micron of PCL length
Globular_lobulesNorm=Globular_lobules./lobule_lengthMicron*100;
Golgi_lobulesNorm=Golgi_lobules./lobule_lengthMicron*100;
MLI2_lobulesNorm=MLI2_lobules./lobule_lengthMicron*100;
Lugaro_lobulesNorm=Lugaro_lobules./lobule_lengthMicron*100;



figure;
b=bar(nanmean(CC_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=color6;
hold on;
for SliceFolder=1:numel(HCR)
    if isempty(strfind(HCR(SliceFolder).slice,'b3'))
     markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
scatter([1:size(CC_lobulesNorm,1)] ,CC_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end
ylim([0 3])
title('CC norm')
box off
set(gcf,'color','none');
set(gca,'color','none');

figure;
b=bar(nanmean(Globular_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=color7;
hold on;
for SliceFolder=1:numel(HCR)
    
     if isempty(strfind(HCR(SliceFolder).slice,'b3'))
     markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
scatter([1:size(Globular_lobulesNorm,1)] ,Globular_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end
ylim([0 3])
title('Globular norm')
box off
set(gcf,'color','none');
set(gca,'color','none');


figure;
b=bar(nanmean(MLI2_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=[color2];
b.FaceAlpha=0.6
hold on;
for SliceFolder=1:numel(HCR)
     if isempty(strfind(HCR(SliceFolder).slice,'b3'))
     markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
scatter([1:size(MLI2_lobulesNorm,1)] ,MLI2_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end

title('MLI2 norm')
box off
set(gcf,'color','none');
set(gca,'color','none');


figure;
b=bar(nanmean(Golgi_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=[177 179 182]./255;
hold on;
for SliceFolder=1:numel(HCR)
     if isempty(strfind(HCR(SliceFolder).slice,'b3'))
     markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
scatter([1:size(Golgi_lobulesNorm,1)] ,Golgi_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end
ylim([0 3])
title('Golgi norm')
box off
set(gcf,'color','none');
set(gca,'color','none');


figure;
b=bar(nanmean(Lugaro_lobulesNorm,2))
b.EdgeColor='none'
b.FaceColor=[127 63 152]./255; %lugaro color
b.FaceAlpha=0.6;
hold on;
for SliceFolder=1:numel(HCR)
    if isempty(strfind(HCR(SliceFolder).slice,'b3'))
     markerColor= [0.5 0.5 0.5];
    else
        markerColor=[0 0 0];
    end
scatter([1:size(Lugaro_lobulesNorm,1)] ,Lugaro_lobulesNorm(:,SliceFolder),'filled','MarkerFaceColor', markerColor)
end
ylim([0 3])
title('Lugaro norm')
box off
set(gcf,'color','none');
set(gca,'color','none');


%% plot cell locations in all slices
Golgicolor=[177 179 182]./255;
Lugarocolor=[127 63 152]./255;
MLI2color=color2;
CCscolor=color6;
Globularcolor=color7;

alphaPCL=0.1;
alphaML=0.7;
mksize=2.5;
lobuleIdx=[1:8]
fig1=figure(1);
fig1.Renderer='Painters';
set(gcf,'position',[680 63 745 915])

for sliceNumber=1:numel(HCR)
    
subaxis(numel(HCR),5,1+((sliceNumber-1)*5), 'SpacingVert',0,'spacingHoriz',0)
 %CCs
Cells=HCR(sliceNumber).CC(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;

% plot shit
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
temp=[1 cumsum(HCR(sliceNumber).CCsperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',CCscolor,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 
    

%Globular
    
subaxis(numel(HCR),5,2+((sliceNumber-1)*5), 'SpacingVert',0,'spacingHoriz',0)
Cells=HCR(sliceNumber).Globular(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;


% plot shit
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
temp=[1 cumsum(HCR(sliceNumber).Globularperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',Globularcolor,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 

%Golgi  
subaxis(numel(HCR),5,3+((sliceNumber-1)*5), 'SpacingVert',0,'spacingHoriz',0)
Cells=HCR(sliceNumber).Golgi(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;


% plot shit
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
temp=[1 cumsum(HCR(sliceNumber).Golgiperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',Golgicolor,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 
 

%Lugaro
subaxis(numel(HCR),5,4+((sliceNumber-1)*5), 'SpacingVert',0,'spacingHoriz',0)
Cells=HCR(sliceNumber).Lugaro(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;

%plot shit
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);

temp=[1 cumsum(HCR(sliceNumber).Lugaroperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',Lugarocolor,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 

%MLI2
subaxis(numel(HCR),5,5+((sliceNumber-1)*5), 'SpacingVert',0,'spacingHoriz',0)
if ~isempty(HCR(sliceNumber).MLI2)
    
Cells=HCR(sliceNumber).MLI2(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;

%plot shit
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
temp=[1 cumsum(HCR(sliceNumber).MLI2perlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',MLI2color,'MarkerEdgeColor','none','MarkerFaceAlpha',1);
end
else 
end
xlim([900 12700])

box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 
    
    
end





%% Make same plot as above but collapse all slices into one.


alphaPCL=0.1;
alphaML=0.7;
mksize=2.5;
mkalpha=0.6;
lobuleIdx=[1:8]
fig2=figure(2);
fig2.Renderer='Painters';
set(gcf,'position',[654 68 248 879])

for sliceNumber=1:numel(HCR)
    
subaxis(5,1,1, 'SpacingVert',0,'spacingHoriz',0)
 %CCs
Cells=HCR(sliceNumber).CC(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;

% plot shit
if sliceNumber==2
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
else 
end
temp=[1 cumsum(HCR(sliceNumber).CCsperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',CCscolor,'MarkerEdgeColor','none','MarkerFaceAlpha',mkalpha);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 
    

%Globular
    
subaxis(5,1,2, 'SpacingVert',0,'spacingHoriz',0)
Cells=HCR(sliceNumber).Globular(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;

hold on
% plot shit
if sliceNumber==2
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
else
end
temp=[1 cumsum(HCR(sliceNumber).Globularperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',Globularcolor,'MarkerEdgeColor','none','MarkerFaceAlpha',mkalpha);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 

%Golgi  
subaxis(5,1,3, 'SpacingVert',0,'spacingHoriz',0)
Cells=HCR(sliceNumber).Golgi(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;

hold on
% plot shit
if sliceNumber==2
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
else
end
temp=[1 cumsum(HCR(sliceNumber).Golgiperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',Golgicolor,'MarkerEdgeColor','none','MarkerFaceAlpha',mkalpha);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 
 

%Lugaro
subaxis(5,1,4, 'SpacingVert',0,'spacingHoriz',0)
Cells=HCR(sliceNumber).Lugaro(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;
hold on
%plot shit
if sliceNumber==2
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
end
temp=[1 cumsum(HCR(sliceNumber).Lugaroperlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',Lugarocolor,'MarkerEdgeColor','none','MarkerFaceAlpha',mkalpha);
end
xlim([900 12700])
box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 

%MLI2
subaxis(5,1,5, 'SpacingVert',0,'spacingHoriz',0)
if ~isempty(HCR(sliceNumber).MLI2)
    
Cells=HCR(sliceNumber).MLI2(:,[1:2]);
PC=HCR(sliceNumber).PC;
ML=HCR(sliceNumber).ML;

%plot shit
if sliceNumber==2
plot(PC(:,1),PC(:,2),'color', [0 0 0 alphaPCL]);
hold on
plot(ML(:,1),ML(:,2),'color', [0 0 0 alphaML]);
else 
end
temp=[1 cumsum(HCR(sliceNumber).MLI2perlobule)+1];
for l=2:numel(temp)
    
    scatter(Cells([temp(l-1):temp(l)],1),Cells([temp(l-1):temp(l)],2),mksize,'MarkerFaceColor',MLI2color,'MarkerEdgeColor','none','MarkerFaceAlpha',mkalpha);
end
else 
end
xlim([900 12700])

box off
set (gca,'Ydir','reverse')
a=600
ylim([400 9000])
% axis square
axis off 
    
    
end

